转录组定量可以用替换featureCounts代替HTSeq-count

在转录组定量分析时，如果采用的是alignment-based转录组定量策略，那么一般会使用的是HISAT2、STAR或者TopHat等比对软件。

接着则是对转录组进行定量，如果是基于基因水平的定量，我之前一般是采用HTSeq-count工具来获取每个基因上的count数。所谓count数，个人简单的理解为根据不同比对情况，将reads分配到各个基因上。HTSeq-count对于多重比对的reads则是采取舍弃策略。当HTSeq-count选择默认参数（-m 默认模式），那么reads是以下图所示的union的情况进行分配的

除了HTSeq-count工具外，其实也可以使用bedtools工具的multicov进行简单的基因水平定量。其需要一个所有基因的位置信息的bed文件，然后计算比对结果bam文件中的reads出现在基因interval上的个数，功能比较简单（说白了就是基于reads位置信息），适用性和准确性一般来说是没有HTSeq-count好。

以上两款在基因水平定量的软件有一个共同点，就是慢！如果样本数上升到几十及上百个样时，其所消耗的时间是以天记的，这时则可以考虑用featureCounts这款软件了，缺点还不知道，但是其优点已经听到过N遍了，就是快。。。

软件下载及安装

下载肯定去官网下载http://subread.sourceforge.net/

featureCounts: a ultrafast and accurate read summarization program

现在featureCounts已经整合在Subread里面了，粗略看了下简介和pdf文档，发现Subread是个功能很全面的软件，而且还有相对应的R包Rsubread，有机会再去瞧瞧。Subread有二进制版本，那么直接下载解压即可使用了

软件的使用

软件的作者认为其软件的优点在于（我就复制黏贴了）：

• It carries out precise and accurate read assignments by taking care of indels, junctions and structural variants in the reads

• It takes only half a minute to summarize 20 million reads（真是快。。。）

• It supports GTF and SAF format annotation

• It supports strand-specific read counting

• It can count reads at feature (eg. exon) or meta-feature (eg. gene) level

• Highly flexible in counting multi-mapping and multi-overlapping reads. Such reads can be excluded, fully counted or fractionally counted（这点跟HTSeq-count不一样了，其对于多重比对的reads并不是只采用全部丢弃的策略，按照其说法是更加灵活的对待）

• It gives users full control on the summarization of paired-end reads, including allowing them to check if both ends are mapped and/or if the fragment length falls within the specified range（可让使用者更加个性化的使用）

• Reduce ambuiguity in assigning read pairs by searching features that overlap with both reads from the pair

• It allows users to specify whether chimeric fragments should be counted（考虑的有点周到）

• Automatically detect input format (SAM or BAM)

• Automatically sort paired-end reads. Users can provide either location-sorted or namesorted bams files to featureCounts. Read sorting is implemented on the fly and it only incurs minimal time cost

featureCounts的参数有点多，可用-h先查看下

1. ~/biosoft/subread/subread-1.6.0-Linux-x86_64/bin/featureCounts -h

按照官网教程简单使用下：

1. ~/biosoft/subread/subread-1.6.0-Linux-x86_64/bin/featureCounts -T 6 -p -t exon -g gene_id -a ~/annotation/mm10/gencode.vM13.annotation.gtf -o SRR3589959_featureCounts222.txt SRR3589959.bam

主要的参数：

-a 输入GTF/GFF基因组注释文件

-p 这个参数是针对paired-end数据

-F 指定-a注释文件的格式，默认是GTF

-g 从注释文件中提取Meta-features信息用于read count，默认是gene_id

-t 跟-g一样的意思，其是默认将exon作为一个feature

-o 输出文件

-T 多线程数

其他参数介绍只能看文档了，不常用的话也是记不住的，要用 49 30754 49 15264 0 0 1126 0 0:00:27 0:00:13 0:00:14 2930再去翻就行

运行中和运行后有两张图可以看看，主要讲了其运行中的一些信息，如下：

featurecount2

从Multimapping reads和Multi-overlapping reads均为not counted可看出我的运行参数是不考虑multimapping的情况的，并且如果一个read在多个基因上有overlaps(multi-overlapping)，则软件也是不计入count数的（只对RNA-SEQ数据，ChIP-seq数据则是例外）

Min overlapping bases为1则说明--minOverlap参数默认是1，因而只考虑全部overlap的reads

Successfully assigned fragments : 22730409 (59.9%)则说明只有59.9%的paired reads定量到了基因上，如果想具体看看那些没有分配上的reads是出于原因，则可看输出中的summary文件

featurecount3

上图的相关说明可以看软件的pdf文档，其实也蛮好理解的

Running time : 0.33 minutes只能说真是快。。。

SRR3589959_featureCounts222.txt文件包含了featureCount结果的信息，如果想了解每个基因上的count数，则只需要提取出第1列和第7列的信息

1. cut -f 1,7 SRR3589959_featureCounts.txt |grep -v '^#' >feacturCounts.txt

然后统计下featureCounts结果总count数和HTSeq-count结果的总count数（这个之前转录组学习时已经跑过了）

1. perl -alne '`$sum += $`F[1]; END {print $sum}' feacturCounts.txt

2. #22730409

3. grep -v '^_' HTSeq-count.txt |perl -alne '`$sum += $`F[1]; END {print $sum}'

4. #21514247

从结果上看，两者数量相差还是蛮小的，看样子两者软件的默认参数的对待reads的方式还是差不多的（主要featureCounts的默认参数对待Multimapping reads是相同的，因为从图3可看出，MultiMapping的reads还是蛮多的）

最后用R作个散点图看看两者同个基因的reads数目大小的分布（取log2），并计算个相关系数

1. library(ggplot2)

2. data <- read.table(file = "count.txt", sep = "\t", header = T, row.names = 1, stringsAsFactors = F)

3. data_matrix <- log2(data + 1)

4. model <- lm(featurecounts ~ htseqcount, data_matrix)

5. r2 <- summary(model)$r.squared

6. r2 <- paste0("r^2=", r2)

7. p <- ggplot(data_matrix, aes(x = featurecounts, y = htseqcount)) +

8.  geom_point(colour = "grey60", size = 1) +

9.  theme_light() +

10.  stat_smooth(method = lm, se = FALSE) +

11.  annotate("text", label = r2, x = 3.5, y = 1)

12. p

从下图中可看出，两者的数据是呈正相关，两者绝大部分的counts数是非常接近的，在低表达量的那部分数据中HTSeq-count的值偏小于featureCounts。因而至少这个数据来说，我觉得featureCounts代替HTSeq-count进行表达量定量是完全没问题的。

featureCount_htseq-count